Time-step targetting methods for real-time dynamics using DMRG 



Adrian E. Feiguin and Steven R. White 

Department of Physics and Astronomy 
University of California, Irvine, CA 92697 
(Dated: February 2, 2008) 

We present a time-step targetting scheme to simulate real-time dynamics efficiently using the 
density matrix renormalization group (DMRG). The algorithm works on ladders and systems with 
interactions beyond nearest neighbors, in contrast to existing Suzuki- Trotter based approaches. 

PACS numbers: 71.27.+a, 71.10.Pm, 72.15.Qm, 73.63.Kv 



Over the last ten years the density matrix renormal- 
ization group (DMRG) Q has proven to be remarkably 
effective at calculating static, ground state properties 
of one-dimensional strongly correlated systems. During 
this period there has also been substantial progress made 
in calculating frequency dependent spectral functions 
However, the most significant progress in extending 
DMRG since its invention has occurred in the last year 
or two. Through a convergence of quantum informa- 
tion and DMRG ideas and techniques, a number of new 
approaches are being developed. The first of these are 
highly efficient and accurate methods for real-time evolu- 
tion, allowing both the calculation of spectral functions 
via Fourier transforming, and also novel time develop- 
ment studies of systems out of equilibrium. 

The key real-time methods thus far developed 0, 0, Q 
rely on the Suzuki- Trotter (S-T) break-up of the evolu- 
tion operator. This approach has a number of impor- 
tant advantages: it is surprisingly simple and easy to 
implement in an existing ground state DMRG program; 
the time evolution is very stable and the only source of 
non-unitarity is the truncation error; and the number of 
density matrix eigenstates needed for a given truncation 
error is minimal. It also has two notable weaknesses: 
it has an error proportional to the time step r squared, 
and, more importantly, it is limited to systems with near- 
est neighbor interactions on a single chain. As we show 
here, the accuracy can be improved using higher order 
expansions. The nearest-neighbor/single chain limitation 
is more problematic. In the case of narrow ladders with 
nearest-neighbor interactions, one can avoid the problem 
by lumping all sites in a rung into a single supersite. 
Unfortunately, this approach becomes very inefficient for 
wider ladders, and is not applicable to general long-range 
interaction terms. 

In this paper we propose a new time evolution scheme 
which produces a basis which targets the states needed 
to represent one time step. Once this basis is complete 
enough, the time step is taken and the algorithm pro- 
ceeds to the next time step. This targetting is intermedi- 
ate to previous approaches: the Trotter methods target 
precisely one instant in time at any DMRG step, while 
Luo, Xiang, and Wang's approach|j| targetted the en- 
tire range of time to be studied. Targetting a wider 



range of time requires more density matrix eigenstates 
be kept, slowing the calculation. By targetting only a 
small interval of time, our approach is nearly as efficient 
as the Trotter methods. In exchange for the small loss 
of efficiency, we gain the ability to treat longer range in- 
teractions, ladder systems, and narrow two-dimensional 
strips. In addition, the accuracy is much improved over 
the lowest order Trotter method. 

We want to find the solution to the time-dependent 
Schrodinger equation: 

ij t m)) = {H{t)-E )m)), a) 

where the ground state energy Eg is introduced to re- 
duce the amplitude of the oscillations by making the di- 
agonal elements of H smaller 0. We use a time de- 
pendent Hamiltonian to include the case where a time- 
dependent perturbation V(t) is added to the time inde- 
pendent Hamiltonian Hq. The initial state \ip(t = 0)) 
is typically the ground state, or the ground state acted 
upon by an operator, but other possibilities are also in- 
teresting. 

When the wavefunction of the system evolves in time, 
its density matrix samples a region of the Hilbert space 
that changes continuously. The DMRG basis is built to 
represent the states that are put into the density matrix 

P = ^2wt\^t)(i>t\, 
t 

where the target states \ipt) are weighted with a factor 
Wt, with J2t w t — 1- Typically, some sweeps are needed 
to build self-consistency between the target states and 
the basis produced by the density matrix. A notable 
exception to this need for self-consistency are the Trotter- 
based time evolution methods: the bond time-evolution 
operator is represented exactly in the current basis, and 
so the pretruncation density matrix is exact. Thus, the 
truncation error is an exact measure of the error in the 
basis produced at that step. In our time-step targetting 
approach, this ideal behavior is lost, and a sweep or two 
is needed to produce a good basis for the time-step. 

How do we produce a density matrix representing the 
wavefunction over an interval of time? Luo, Xiang and 
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Wang (see also Ref. @) suggested targetting the 
wavefunction at a sequence of times spanning the inter- 
val, ijj{t = 0), ip(t — t), ip(t = 2r), . . ., ip{t — tit), simul- 
taneously. We argue that this choice is very close to ideal. 
Suppose that our basis includes ip{kr) and ip({k + l)r). 
Then the basis includes any linear combination of these 
states, so that one could imagine using an interpolation 
formula to determine coefficients a and b to approximate 
the wavefunction at any time between kr and (k + l)r as 
ip(t) w a^{kr) + bip((k + l)r). This suggests that the er- 
ror in the basis is at worst r 2 . If the basis includes more 
than two time points, one could imagine using higher or- 
der interpolations, e.g. splines, putting a tighter bound 
on the error in the basis. The key point is that we do not 
actually perform these interpolations; the basis is auto- 
matically good enough to allow whatever interpolation is 
most accurate given the set of time points. This suggests 
that the error in the basis varies as r™ +1 . 

If r is small enough and n big enough, and enough 
self-consistency sweeps are made, the error in the basis is 
given by the truncation error. This is the ideal situation 
for a DMRG calculation. Since this truncation error is of- 
ten miniscule we shall say that an approximate algorithm 
is "quasiexact" when the error is strictly controlled by the 
DMRG truncation error e (with some properties propor- 
tional to e and other to e 1 / 2 ). For example, the infinite 
system method applied to a finite system is not quasiex- 
act, even though the error goes to zero as the number of 
states kept m increases. If enough sweeps are taken, and 
absent any "sticking" problems with metastable ground 
states, the finite system ground state DMRG method is 
quasiexact. Non quasiexact algorithms seem to be the 
source of most DMRG "mistakes" . The procedure below 
is nearly quasiexact: it has a small separate time step 
error, perhaps of order r 4 , in addition to the truncation 
error. 

Our procedure consists of taking a tentative time step 
at each DMRG step, the purpose of which is to generate a 
good basis. The standard fourth order Runge-Kutta (R- 
K) algorithm is very convenient for this purpose. This is 
defined in terms of a set of four vectors: 

I**) = T H(t)\m), 

|*2> = rH{t + r/2) [|^(t)> + l/2|fci>] , 
|fc 3 ) = TH(t + T/2)[\i;(t)) + l/2\k 2 )}, 
\k 4 ) = TH(t + T)[\ij(t)) + \h)}, (2) 

where H{t) — H(t) — E$. The state at time t + t is given 
by 

I-0C* + r)) « \ \\ki) + 2\k 2 ) + 2\k 3 ) + |fc 4 }] + 0(r 5 ). (3) 
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We choose to target the state at times t, t+r/3, t + 2r/3 
and t + r. The R-K vectors have been chosen to minimize 
the error in \ip(t + r)), but they can also be used to 
generate \tfj) at other times. The states at times t + r/3 
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FIG. 1: Error E(t = 8) for the Haldane chain (L = 32), 
according to Eq. JSJ, as a function of the number of states 
kept m. We show results of simulations using the Runge- 
Kutta algorithm, for different time steps and distributions of 
weights. 



and t + 2t/3 can be approximated, with an error 0{t ), 
as 



m+r/3)) a \m) 
i 



+ -^[3I|fc 1 ) + I4|fc 2 ) + 14|fc 3 )-5|fc4)], 
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m+2r/3)) « + 

+ i [16|*i) + 20|fc 2 ) + 20|fc 3 ) - 2|fc 4 )P) 

In practice we proceed as follows: each half-sweep cor- 
responds to one time step. At each step of the half-sweep, 
we calculate the R-K vectors ©, but without advancing 
in time. The density matrix is then obtained by using 
the formula Q with the target states \tp{t)), |-0(t + r/3)), 
\ip(t + 2r/3)), and \tjj(t + t)). Advancing in time is done 
on the last step of a half-sweep. However, we may choose 
to advance in time only every other half-sweep, or only 
after several half-sweeps, in order to make sure the basis 
adequately represents the time-step. Our tests show that 
one half-sweep is adequate and most efficient for the sys- 
tems studied here. The method used to advance in time 
in the last step need not be the R-K method used in the 
previous tentative steps. In fact, the computation time 
involved in the last step of a sweep is typically miniscule, 
so a more accurate procedure is warranted. A simple 
way which keeps the time- integration errors much smaller 
than the basis errors is by performing 10 R-K iterations 
with step r/10. We usually use this method. Alterna- 
tively, one can evolve using the exponential of the the 
Hamiltonian in the Lanczos tridiagonal representation, 
which is exactly unitary. However, the truncation to a 
finite number of density matrix eigenstates introduces 
nonunitarity anyway, so the Lanczos procedure has no 
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FIG. 2: Same as in Fig. using 1st order Suzuki- Trotter 
break-up (gray symbols), 4th order Suzuki- Trotter (empty 
symbols), and 4th order Runge-Kutta (filled symbols). 



FIG. 3: Error E(t = 8) for the Haldane chain for different 
time steps r: a) 1st, 2nd, and 4th order Suzuki- Trotter break- 
ups and m = 160; b) Runge-Kutta and m — 100. 



special advantage. In practice, we find comparable over- 
all accuracy in the two methods 0. 

To test the method we first studied the S = 1 Heisen- 
berg chain. Since it is a single chain system, the Suzuki- 
Trotter methods are also applicable. In addition to our 
new method, we have used both the traditional 1st or- 
der Suzuki- Trotter decomposition Q and the 4th order 
Forest- Ruth break-up [l0(. In order to compare the re- 
sults, we calculated the error as 



E(t) = 



\ 



Exact 



(x,t)Y 



(5) 



where S z Exact is obtained using 4th order Suzuki- Trotter 
with m = 200 and r = 0.02, which keeps the truncation 
error under 10~ 12 . 

The target states (@J can be weighted equally, or un- 
equally. We have performed several test runs with differ- 
ent distributions of weights. In Fig. ^ we show the error 
(J5J at time t = 8 as function of the number of states kept 
m, for various weightings. The best weighting we have 
found is wi — w 3 — 1/3, w 2 — w 3 — 1/6. The calcu- 
lations described below, unless otherwise noted, use this 
choice of weights. 

In Fig. [2] we compare results by using our method 
and Suzuki- Trotter evolution. The Suzuki- Trotter simu- 
lations converge when the error reaches a plateau and re- 
mains constant with increasing number of states m. This 
occurs generally for a relatively small m, after which the 
accuracy of the simulation is completely controlled by 
the Trotter error, and not by the truncation error. In 
Fig. we verify that the quantity E is proportional 
to r, t 2 , and t 4 for the three Suzuki- Trotter break-ups 
considered [lo|. In the R-K simulations, convergence is 
slower with the number of states m because we need the 



basis to be optimized for 4 states at slightly different 
times. The accuracy improves steadily with the size of 
the basis, and also with the size of the time-step r, as it 
can be seen in Fig. |3Jd, although the method breaks down 
for time-steps larger than t ~ 0.25. These results may 
look counter intuitive, since the R-K error is expected 
to be proportional to t 4 . The reason for this behaviour 
is that smaller time-steps require more iterations, with 
a consequent accumulation of error due to the trunca- 
tion. Therefore, unlike the S-T case, the simulation is 
now dominated by the truncation error, which can be 
reduced by increasing the size of the DMRG basis. 

For typical accuracies on a ID chain, the R-K algo- 
rithm is numerically costlier than the S-T counterparts. 
Measuring the the CPU time required to reach a time t 
in the simulation we find, for instance, that in order to 
obtain an error of the order of 10 -3 at t = 8 we could 
use 1st order S-T with r = 0.016 and m = 40 (CPU time 
for 1000 half-sweeps: 34 minutes), 4th order S-T with 
t = 0.25 and m = 40 (CPU time for 224 half-sweeps: 7 
minutes), or R-K with r = 0.10 and m = 140 (CPU time: 
104 minutes). The R-K method requires fewer sweeps to 
reach a specified time. However, it also requires more 
states to be kept, leading to a substantially lower com- 
putation time. 

In Fig. ^ we show how the number of states m re- 
quired to keep a fixed, very small truncation error of 
10~ 8 grows with time. This rapid growth in m for a 
fixed accuracy is not surprising. At t = 0, an operator 
is applied to the ground state, creating \ip(Q)). For small 
t, \tp(t)) is still closely related to the ground state, and 
so requires a comparable number of states to represent 
it. For larger t, |^>(t)) becomes more complicated as each 
excited eigenstate evolves with a different frequency and 
becomes independent of the others. 
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FIG. 4: Number of states required to keep a truncation error 
of 1CP 8 , as a function of time. The results correspond to a 
R-K simulation of a Haldane chain with L — 32. 



As an application of the R-K method we calculated 
the spin structure factor for a 2 x L Heisenberg ladder 
with spin S — 1/2, A(k,w) = — ^ImG(k, uj), obtained by 
Fourier transforming the time-dependent spin-spin cor- 
relation function 

G(x,i) = (S-(x,t)5 + (0,0)). 

In this case, besides targetting the four states at different 
times (0J , we also need to target the ground state at t = 
0. We have used a weight wo — 1/2 for the ground state, 
and all the other weights equal to 1/8. In Fig. |SJwe show 
the results for L = 32 using a time step r = 0.1 and m = 
256, which kept the truncation error under 1CP 7 for times 
up to t — 30. The result for the spin-gap is A = 0.506, 
which should be compared to the very precise DMRG 
value AExact = 0.50249 in the thermodynamic limit [llj |. 
We also show for comparison the exact diagonalization 
results for the singlet-triplet excitations for L = 12 from 
Ref.^3- A continuum of excitations can be observed 
above the magnon band for k y = 0. It becomes more 
difficult to resolve the band for k y = in the proximity 
of k x — > because the quasiparticle weight tends to zero 
in this limit. This is not the case for k y — n, where the 
band is well defined in the entire range of momenta. 

To summarize, we have presented a new algorithm for 
simulating time evolution of quantum systems. We de- 
scribed how to tune the parameters in order to reach 
accuracies comparable to those obtained by using Susuki- 
Trotter based approaches, and demonstrated its applica- 
tion by calculating the excitation spectrum of the Heisen- 
berg ladder. Unlike methods that rely on Suzuki- Trotter 
break-ups, our algorithm can be applied to systems with 
arbitrary geometry, and interactions beyond first neigh- 
bors. Moreover, it can be easily generalized for studying 



more complex models with strong correlations. 
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FIG. 5: Structure factor A(k,u)) for the Heisenberg ladder 
using 4th order Runge-Kutta, m = 256 states, and time-step 
r = 0.1. The solid line is centered atthe quasiparticle peak. 
The tones of gray are proportional to the quasiparticle weight 
(dashed curve). The symbols are Lanczos results for L = 12 
from Ref.H3]. 
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